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ABSTRACT 

A  combined  theoretical  and  empirical  numerical  model  was 
developed  which  predicts  the  performance  of  continuous 
electrode  coaxial  magnetoplasmadynamic  thrusters  as  a 
function  of  thruster  dimensions,  mass  flow  rate,  and  input 
current.  This  model  was  used  to  predict  the  effects  of 
scaling  on  these  thrusters. 

The  model  predicts  that  for  scaling  factors  down  to  one- 
half,  relations  can  be  found  relating  the  performance  of  one 
thruster  to  another.  The  model  was  used  to  examine  these 
relationships  for  four  different  thruster  configurations 
over  a  broad  range  of  operating  currents.  The  thrusters 
examined  consisted  of  two  geometries  and  their  half  scale 
counterparts.  A  conclusion  from  the  analysis  is  that 
scaling  down  the  size  of  the  thruster  by  50%  can  reduce  the 
total  power  input  by  30%  to  40%  at  comparable  efficiencies. 
However,  this  is  at  the  cost  of  increasing  the  specific 
impulse  by  a  factor  of  two  which  may  render  the  thruster 
inappropriate  for  the  intended  missions. 
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I.  INTRODUCTION 

Equation  1  expresses  the  relationship  between  initial  and 
final  space  vehicle  mass  to  the  propellant  exhaust  velocity 
and  the  velocity  increment  required  for  a  maneuver  in  space 
[1].  This  equation  shows  that  a  rocket  engine's  exhaust 
velocity  must  be  greater  than  the  velocity  increment  of  the 
mission  of  interest  for  a  significant  amount  of  payload  to 
be  delivered.  If  the  exhaust  velocity  is  much  less  than  the 
velocity  increment,  the  ratio  of  spacecraft  masses  becomes 
small  implying  a  small  final  mass  or  payload.  The  velocity 
increment  depends  only  on  flight  path  and  is  independent  of 
the  vehicle.  Therefore,  the  only  way  of  improving  the  mass 
ratio  is  to  increase  the  exhaust  velocity  of  the  rocket 
engine  . 


Mf/Mi  =  exp(-Av/ue) 


(1) 


Chemical  rocket  engines  have  a  maximum  exhaust  velocity 
determined  by  the  amount  of  energy  stored  in  the  chemical 
bonds  of  the  fuel.  This  limits  exhaust  velocities  to  less 
than  5000  meters  per  second.  Orbit  transfer  missions  such 
as  Low  Earth  Orbit  (LEO)  to  Geosynchronous  Orbit  (GEO) 
require  velocity  increments  of  approximately  4200  meters  per 
second.  Application  of  these  values  to  Equation  1  gives  a 
mass  ratio  of  less  than  one-half.  This  implies  that,  at 
best,  less   than  half   of  the  initial  weight  of  a  spacecraft 


can  be  useful  payload.  For  current  space  motors,  the  ratio 
of  initial  mass  to  final  mass  is  about  4:1.  An  alternative 
to  chemical  propulsion  is  electric  propulsion.  Electric 
propulsion  avoids  the  thermodynamic  limitations  of  chemical 
engines  by  employing  an  external  power  source  to  generate 
the  forces  required  to  accelerate  a  propellant.  The  amount 
of  energy  available  to  provide  thrust  is  no  longer 
determined  by  the  nature  of  the  propellants.  This  enables 
greater  amounts  of  energy  to  be  applied  to  the  propellant 
than  possible  with  chemicals.  One  such  device  is  the 
magnetoplasmadynamic  (MPD)  thruster.  This  coaxial  electric 
propulsion  device  uses  the  Lorentz  jxB  force  to  accelerate  a 
plasma  to  velocities  in  excess  of  25,000  meters  per  second. 
This  increased  exhaust  velocity  enables  a  50%  to  200% 
increase  in  the  mass  of  the  delivered  payload.  An  electrical 
power  supply  must  now  be  included  in  the  spacecraft  mass, 
but  the  reduction  in  propellant  mass  enabled  by  the  MPD 
thruster's  high  exhaust  velocities  may  be  suffucient  to 
allow  for  the  mass  of  the  power  supply.  Figure  1  shows  an 
MPD  thruster. 


Figure  1   Continuous  Electrode  MPD  Thruster 


It  can  be  seen  that  the  thruster  is  structurally  a  very 
simple  device.  The  MPD  thruster  has  attracted  considerable 
attention  not  only  for  its  simplicity,  but  also  because  much 
greater  thrust  levels  than  any  other  electric  propulsion 
device  can  be  generated.  This,  combined  with  its  ability  to 
operate  with  virtually  any  propellant,  make  it  a  prime 
candidate  for  supporting  space  missions. 

The  price  to  be  paid  for  eliminating  the  dependence  on 
the  chemical  energy  in  the  fuel  and  oxidizer  is  the 
requirement  of  an  external  power  source.  Power  requirements 
for  even  modest  sized  MPD  thrusters  reach  to  the 
mul t imegawat t  range. 


Space  based  power  supplies  of  this  size  will  not  be 
available  for  many  years.  Two  options  are  available  to 
reduce  the  power  requirements  of  the  thruster.  The  first  is 
to  increase  the  efficiency  of  the  thruster  design.  Since 
these  thrusters  are  now  operating  in  the  50%  efficiency 
regime  [2],  this  technique  promises  a  reduction  in  power 
requirements  of  less  than  a  factor  of  two.  The  second 
option  is  to  develop  smaller  thrusters  which  require  less 
power  than  their  full  size  counterparts.  MPD  thrusters,  by 
virtue  of  their  high  thrust  density,  would  still  be  able  to 
produce  useful  levels  of  thrust  after  significant  scaling. 
This  technique  promises  much  greater  factors  of  power 
reduction  than  the  first  option  because  power  requirements 
are  related  to  size  by  a  term  of  greater  than  linear  scale. 

It  is  the  goal  of  this  thesis  to  develop  a  computer  model 
of  an  MPD  thruster  which  will  predict  specific  impulse  and 
efficiency  as  functions  of  current,  thruster  dimensions,  and 
propellant  mass  flow  rate. 


II.  DERIVATION  OF  THE  MPD  THRUST  MODEL 

Equations   2,  3,  and  4  are  the  equations  of  conservation 
of  mass,  momentum,  and  energy  respectively. 


3p/3t  +  V'(pu)  =  0 


(2) 


p[3u/3t  +  (u'V)u]  =  -Vp  +  (j*B)  +  fv  (3) 


p(3  /at  +  u-V)(cpT  +  u2/2)  =  3p/3t  +  j'E  +  4>£  +  4>v  -  4>r   (4) 

Since  an  MPD  thruster  is  intended  to  be  operated  in  a 
continuous  mode,  these  equations  can  be  simplified  to  the 
steady  state  forms  of  Equations  5,  6,  and  7.  This  is 
justified  because  measurements  of  the  start  up  and  shut  down 
transients  in  MPD  thrusters  are  measured  in  seconds  where  as 
actual  thruster  operation  is  intended  to  last  for  up  to  100 
days  [ 3 ] . 


V-(pu)  =  0 


(5) 


p(TI'V)u  =  -Vp  +  (jxB)  +  fv 


(6) 


p(u-V)(cpT    +    u2/2)    =    j-E    +    4>£    +    <J>V    -    4>r 


(7) 


For  the  derivation  of  an  expression  for  the  thrust  of  an 
MPD  thruster  Equation  6,  conservation  of  momentum,  will  be 
of  paramount  importance.  The  low  density  and  viscosity  of 
the  gaseous  argon,  which  will  be  the  propellant  used,  within 
the  arc  chamber  suggests  that  the  viscous  drag  density  can 
be  assumed  negligible.  Calculations  [4]  show  that  this 
force,  when  integrated  over  a  typical  thruster,  is  only  on 
the  order  of  a  tenth  of  a  Newton  where  as  the  total  thrust 
is  typically  three  orders  of  magnitude  greater  than  this. 
This  simplification  results  in  Equation  8. 


p(u«V)u  =  -Vp  +  (jxB) 


(8) 


It  is  assumed  the  thruster  in  question  is  a  small  section 
of  an  infinitely  long  coaxial  device.  This  assumption 
greatly  simplifies  the  mathematics  by  eliminating  the 
requirement  to  deal  with  electric  and  magnetic  field 
fringing  and  the  exit  plane  of  the  thruster.  In  addition  it 
will  be  assumed  that  the  thruster  consists  of  a  series  of 
coaxial  streamtubes   of  cylindrical   cross  section   and  that 


steady  state  migration  of  propellant  across  the  stream  tube 
boundaries  does  not  occur.  This  permits  the  application  of 
one  dimensional  analysis  to  the  problem  of  modeling  the 
forces  generated  by  the  thruster.  As  will  be  seen  later, 
this  division  of  the  thruster  into  coaxial  streamtubes  is 
also  necessary  for  the  application  of  a  numerical 
integration  program. 

The  above  assumptions  and  simplifications  permit  Equation 
8  to  be  expanded  using  cylindrical  coordinates  to  the  form 
of  Equation  9  . 

p[  (ur  +  u(i)  +  nz)  •  (r3  /3r+$3  /3<j>  +  z3  /3z)]'[ur  +  u^+irz]  = 
-[f3p/3r  +  $3p/3(t.  +  23p/3z]  +  [(jr  +  j(t)  +  jz)(Br  +  B(t)  +  Bz)]       (9) 

From  symmetry,  the  assumption  that  the  only  velocity 
which  can  be  maintained  is  in  the  axial  or  "z"  direction, 
and  the  definition  of  the  vector  dot  product,  Equation  9, 
can  be  simplified  to  Equation  10. 


zpu23uz/3z  =  -f3p/3r  -  z3p/3z  +  zjrB^  -  rj^B^ 


(10) 


By   equating  vector  components,  Equation  10  can  be  split 
into  Equations  11  and  12. 


puz3uz/3z  =  -3p/3z  +  j'j-Ba 


(11) 


0  =-3p/3r  -  jzB$ 


(12) 


After   rearranging,  Equation  11  can  be  put  into  the  form 
of  Equation  13  and  Equation  12  becomes  Equation  14. 


fz  =   uz9uz/3z  +  9p/8z  =  jrB( 


(13) 


fr  =  -9p/3r  =  j  B 
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(14) 


It  is  convention  [1],  [5],  [6]),  [7],  [8],  [9]  to  refer 
to  the  f  component  as  the  blowing  force  density  and  the  fr 
component  as  the  pumping  force  density.  The  thrust  produced 
by  an  MPD  thruster  is,  in  large  part,  the  sum  of  these  two 
components  integrated  over  the  volume  of  the  thruster.  The 
remaining  thrust  is  the  result  of  electrothermal  heating  of 
the  propellant  as  it  flows  through  the  arc  chamber. 
Expressions  for  these  components  of  thrust  will  be  developed 
individually  after  which  the  three  will  be  summed  to  give 
the  total  thrust. 

An  mathematical  model  for  thrust  can  only  be  accomplished 
if  the  current  density  within  the  arc  chamber  can  be 
mathematically  defined.  It  has  been  experimentally 
determined  [6]  that  over  a  wide  range  of  thruster  shapes  and 
sizes,  the  current  density  follows  an  "S"  shaped  path  as 
illustrated  in  Figure  2. 
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Figure  2   Current  Flow  Pattern  Inside  an  MPD  Thruster 

For  a  continuous  electrode  thruster,  the  depth  of  this 
curve  or  amount  of  axial  displacement  is  determined  by  the 
Hall  parameter.  For  typical  operating  regimes,  the  Hall 
parameter  varies  from  effectively  zero  to  a  value  of 
approximately  ten  [5],  [6]  depending  on  the  current  for  a 
fixed  thruster.  By  trial  and  error,  Equations  15  to  18 
were  developed  and  provide  the  current  flow  pattern 
required . 


for  rc<r<(rc+ra)/2 

jr=(0.159/z  ){[J/r]2-[2.22xl(I6(r-rc)J2]2}-5 


(15) 


jz  =  3.54xlO"7J2(r-rc)/z' 


(16) 


For  (rc+ra)/2<r<ra 

jr=(0.159/z')([J/r]2-[2.22xl0-6(ra-r)J2]2}-5 


(17) 


jz  =  3.54xl(T7J2(ra-r)/z 


(18) 


With  these  relations  to  define  the  current  within  the 
thruster,  expressions  can  be  found  for  the  blowing,  pumping, 
and  electrothermal  components  of  thrust. 


BLOWING  THRUST 

The  blowing  component  of  thrust  results  from  the  vector 
product  of  radial  current  density  with  an  azimuthally 
oriented  magnetic  field  as  shown  in  Figure  1.  Equations  15 
and  17  provide  the  radial  current  densities.  The  magnetic 
field  is  found  by  applying  Ampere's  law  for  straight 
conductors  . 

B  =  $B  =  $UoJ/27Tr     Ampere's  Law  (19a) 

The  magnetic  field  which  results  is  given  in  equation  (19b). 


10 


B  =  (y0J/2TTr)(l-z/zf  ) 


(19b) 


As  can  be  seen  from  an  examination  of  this  equation,  the 
magnetic  field  is  a  function  both  of  radial  and  axial 
position  within  the  thruster.  Vector  multiplication  of  these 
two  crossed  vectors  results  in  a  vector  directed  in  the 
positive  "z"  direction  with  a  magnitude  defined  by  Equations 
20  and  21. 

For  rc<r<(rc+ra)/2 

fzl  =  (0.159JUo/2Trrz,){[J/r]2-[2.22xlO-6(r-rc)J2]2}-5(l-z/z,)(20) 

For  (rc+ra)/2<r<ra 
fz2=(0.159JUo/2Trrz,){[J/r]2-[2.22xlO-6(ra-r)J2]2}'5(l-z/z,)(21) 

These  last  two  equation-  give  the  force  density,  in  newtons 
per  cubic  meter,  as  a  function  of  radial  and  axial  position. 
To  find  the  total  blowing  thrust  component,  it  is  necessary 
to  integrate  these  force  densities  over  the  internal  volume 
of  the  arc  chamber. 


For  rc<r< ( rr  +  r~ ) / 2 


z'  f(rc+ra)/2 


Fzl  =      jfzldz2Trrdr 
0    r„ 


(22) 
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For    (rc+ra)/2<r<r. 


z2    = 


f  f 

f  22dz27rrdr 

0    (rc+ra)/2 


(23) 


After  integrating   with  respect   to  the   axial  and  azimuthal 
variables,  Equations  24  and  25  result. 


(rr+r„)/2 


zl  = 


c    *  a 


7.95xl0lt([J/r]2-[2.22xl(T6(r-rr)J2]2}-5Jrdr 


(24) 


FZ2: 


7.95xl0M[J/r]2-[2.22xl0-6(ra-r)J2]2}-5Jrdr 
(rc+ra)/2 


(25) 


The  form  of  this  equation  is  such  as  to  resist  integration 
analytically  or  through  the  use  of  generic  forms  given  in 
integral  tables.  The  total  blowing  force  was  obtained  by 
applying  numerical  integration  techniques  to  Equations  24 
and  25.  This  integration  uses  the  program  in  Appendix  B 
titled  Theoretical  Model  Program. 

PUMPING  COMPONENT 

The  pumping  force  density  given  by  Equation  14  is  the 
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result  of  the  vector  product  of  the  azimuthal  magnetic  field 
and  an  axial  component  of  the  current  in  the  negative  "z" 
direction.  This  cross  product  gives  rise  to  a  force  density 
component  in  the  negative  radial  direction.  This  force 
density  is  balanced  by  a  pressure  gradient  in  the  positive 
radial  direction.  Integration  of  the  pressure  gradient  at 
the  z'  =  0  position  yields  an  expression  for  pressure  as  a 
function  of  radial  position.  Multiplying  pressure  with  a 
differential  area  and  integrating  provides  the  total  pumping 
force . 

Expanding  Eequation  14  by  substituting  in  the  definitions 
of  the   axial  current   density  and   the   azimuthal   magnetic 
field  at  the  z  =  0  location  results  in  Equations  26  and  27. 
For  rc<r<(rc+ra)/2 


fr  =-3p/3r  =  cJ3(r-r  )/z'r 


(26) 


For  (rc+ra)/2<r<ra 

fr  =-3p/3r  =  cJ3(ra-r)/z'r 


(27) 


where    c  =  7.08X1011*  newtons/cm2amperes 3 
Separating  variables   and  integrating  gives  Equations  28  and 
29  which   collectively  represent  the  equation  for  the  radial 
pressure  gradient  as  a  function  of  radial  position. 


For  rr<r<(rr+r„)/2 


p  =  p0+c( r-r.lnr ) J3/z 


3  /„  t 


(28) 
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For  (rc+ra)/2<r<ra 

p  =  p0+c(ralnr-r)J3/z'  (29) 

Integrating  Equations  28  and  29,  with  respect  to  a 
differential  ring  area,  and  collecting  terms  yields  equation 
30. 


Fpumping  =  (2TrcJ3/z'){((ra  +  rc)2)3/3-[((ra  +  rc)/2)2 
ln((ra+rc)/2)-((ra+rc)/2)2/4]-[rc3/3- 
rc((rc2/2)lnr_-rc2/A)]  +  [(ra3/2)lnr;,- 


ra3/4-raV3]-[ra(((ra  +  rr)/2)2/2)ln((r 


+rc)/2)-ra((ra+rc)/2)2/4-((r  +rr)/2)3/3]}   (30) 


While  this  expression  is  long,  it  is  algebraically  simple 
and  is  easily  evaluated  numerically.  Like  the  blowing  thrust 
term,  Equation  30  has  been  included  in  the  Appendix  B 
program. 

ELECTROTHERMAL  THRUST 

Due  to  ohmic  heating  the  propellant  in  an  MPD  arc  chamber 
increases    in    temperature    from    an    initial  288  K   to 
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temperatures  exceeding  20,000  K  at  some  locations  within  the 
thruster.  Cory  [4]  derived  an  expression  for  the 
electrothermal  thrust  based  on  the  assumptions  of  negligible 
viscous  forces,  uniform  pressure  on  the  chamber  walls,  and 
no  thruster  pressure  present  on  the  outside  surfaces  of  the 
arc  chamber.  His  result  for  the  pressure  is  presented  here 
as  Equation  31 . 


3  T  1.5  A-3 


p  =  3.02xlO"dJ1'5m 


(31) 


Multiplying  this  expression  by  the  exit  plane  area  of  the 
thruster  provides  an  expression  for  the  electrothermal 
component  of  thrust.  This  expression  is  Equation  32. 


Fe  =  3.02xl03J1-5m3[TT(ra2-rc2)] 


(32) 


TOTAL  THRUST 

As  stated  before,  the  total  thrust  is  the  sum  of  the 
three  thrust  components.  A  Fortran  program  titled 
Theoretical  Model  Program  given  in  Appendix  B  performs  the 
summation.  The  results  of  calculations  using  this  model  are 
discussed  in  section  III. 
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THRUSTER  VOLTAGE 

An  eraperical  expression  for  thruster  voltage  was 
assembled  using  all  available  information  on  coaxial  MPD 
thrusters  tested  to  date,  [2],  (3),  [4],  [5],[6],  [8],  (10), 
and  [11].  First,  an  equation  for  voltage  as  a  function  of 
thruster  length  was  developed.  Second,  the  influence  of 
thruster  diameter  was  included.  Finally,  the  influence  of 
mass  flow  rate  was  incorporated  into  the  equation.  This 
voltage  model  was  tested  using  thruster  geometries  not  used 
in  its  formulation.  The  model  predicted  interpolated 
voltages  to  within  three  percent  and  extrapolated  voltages 
within  seven  percent.  Both  of  these  errors  are  within  the 
plus-or-minus  10  to  30  percent  error  bars  typical  of 
laboratory  testing  associated  with  MPD  thruster  operation. 
Therefore,  the  voltage  model  was  judged  to  be  acceptable. 
This  model  is  represented  by  Equation  33. 

V=( (55.4-26. 31* (l-ra/5.1) -23.4* (l-Z/21. 6) )- 
((39.2-65.8*(l-ra/5.1)-27.7*(l-Z/21.6))*lE-4* 

J+((35.7-54.5*(l-ra/5.1)+3.88*(l-Z/21.6))*lE-8)*J2)* 
(178.0-4  7.3*1000.0*W+4.3*(1000.0*W)2)/49.o  (33) 

In  this  equation,  z  represents  the  thruster  length  and  W  the 
propellant  mass  flow  rate. 
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SUMMARY 

From  the  input  parameters  of  thruster  length  (z),  anode 
radius  (r0),  cathode  radius  (r  ),  current  (J),  and  mass  flow 

3.  C 

rate  (W) ,  the  model  predicts  thrust  and  voltage.  These 
predictions  are  used  in  Equations  34  and  35  to  calculate 
specific  impulse  (Isp)  and  efficiency  (N)  respectively. 

Isp  =  T/(W*9.8)  (34) 

N  =  (0.5*W*(ISP*9.8)2/(V*J)  (35) 

The  model  has  been  written  as  a  Fortran  program  titled 
Theoretical  Model  Program.  A  copy  is  included  in  this 
thesis  as  Appendix  B. 
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III.  APPLICATION  AND  RESULTS  OF  THE  THEORETICAL  MODEL 

The  theoretical  computer  model  was  tested  over  a  wide 
range  of  input  currents  of  known  thrusters.  Graphs  of 
efficiency  versus  specific  impulse  matched  the  experimental 
results  [6]  to  within  plus  or  minus  ten  percent.  Similar 
agreement  was  obtained  for  graphs  of  efficiency  versus  input 
power.  The  plus  or  minus  ten  percent  criteria  was  deemed 
acceptable  because  it  placed  the  model's  plots  within  the 
error  bar  of  the  experimental  data.  This  test  verified  the 
model's  accuracy. 

The  model  was  used  to  predict  the' specific  impulse  and 
efficiency  of  four  different  thrusters  over  a  range  of 
current  levels.  The  thrusters  were  chosen  with  geometries 
which  grouped  them  into  two  pairs  of  full  and  half  sized 
devices.  All  thrusters  were  assumed  to  be  using  argon  as 
propellant.  A  coded  designator  was  created  for  each  device 
to  make  identification  easier.  This  designator  is  of  the 
form  A/B/C/D.  A  is  the  thruster  length.  B  is  the  anode 
radius.  C  is  the  cathode  diameter.  D  is  the  mass  flow  rate 
of  the  argon  propellant  in  kilograms  per  second.  All 
dimensions  are  in  centimeters.  Thus  a  20/5/1/. 006  thruster 
is  20  centimeters  long  with  a  5  centimeter  inside  anode 
radius  and  a  1  centimeter  diameter  cathode.  The  mass  flow 
rate  is  6  grams  per  second.    Thruster  sizes,  aspect  ratios, 
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mass  flow  rates,  and  applied  currents  are  all  consistent 
with  those  regimes  of  operation  typical  of  MPD  devices. 
Very  small  devices  operating  at  power  levels  or  mass  flow 
rates  where  previous  theoretical  or  experimental  evidence 
indicated  operation  would  be  unstable  were  not  examined. 

The  four  thrusters  examined  were  20/5/1/. 006,  10/5/1/. 006 
and  their  half  scale  analogs  10/2 . 5/ . 5/ . 0015 ,  5/2 . 5/ . 5/ . 0015 
respectively.  The  10/5/1/. 006  configuration  was  chosen 
because  it  closely  models  the  Princeton  Benchmark  thruster 
for  which  significant  theoretical  and  experimental  results 
exist.  The  20/5/1/. 006  thruster  was  selected  because  recent 
research  indicates  thrusters  with  extended  electrodes  are 
more  efficient  than  shorter  electrode  thrusters  like  the 
Benchmark.  The  mass  flow  rate  for  the  half  scale  thrusters 
was  reduced  by  the  square  of  the  scale  factor  to  reflect  the 
reduction  in  exhaust  area.  This  maintains  a  constant  atomic 
particle  density.  Failure  to  do  so  would  place  the  scaled 
thrusters  in  a  completely  different  operating  regime. 
Parameters  for  each  of  these  devices  became  the  input  to  the 
computer  model  with  a  range  of  currents  from  1,000  to  25,000 
amperes.  The  model's  outputs  were  plotted  and  are  presented 
in  Figures  3  through  7  inclusive. 

Figures   3   through   5   show   the   relationship   between 
specific  impulse   and  efficiency.    A  comparison  between  the 
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graphs  for  10/5/1/. 006  and  20/5/1/. 006,  Figure  3,  discloses 
that,  as  predicted,  the  longer  configuration  is  more 
efficient  at  a  given  specific  impulse.  The  same  result  is 
observed  in  Figure  4  for  the  half  scale  thrusters.  This  can 
be  explained  by  examining  the  change  in  current  density  and 
current  patterns  as  thruster  lengths  increase.  For  a  given 
current,  a  longer  device  has  a  lower  average  current 
density.  This  results  in  lower  collision  frequencies  and, 
therefore,  smaller  Hall  currents.  Hall  currents,  in  MPD 
thrusters,  travel  in  the  axial  direction  and  are  the  result 
of  the  combined  effects  of  differing  collision  rates  and 
mobilities  between  electrons  and  ions.  It  is  the  Hall 
current  which  creates  the  pumping  component  of  thrust. 
Equation  18  is  an  emperically  derived  expression  for  the 
Hall  current.  Because  this  component  involves  propellant 
being  accelerated  through  radia  pressure  gradients  and  a  90 
degree  turn  to  exit  the  thruster,  it  is  less  efficient  than 
the  blowing  component  which  is  axially  oriented  throughout 
the  acceleration  process.  This  result  could  have  been 
predicted  by  examining  the  relative  magnitudes  of  the  forces 
from  Equations  24,  25,  and  29.  Qualitatively,  both  the 
blowing  and  pumping  thrust  components  are  functions  of  total 
current  raised  to  the  third  power.  However,  the 
multiplicative  constant  for  the  blowing  force  appears  to  be 
larger  than  any  value  expected  from  the  corresponding 
constant  in  the  pumping  equation.  One  would  expect  this 
phenomena  to   continue  until   viscous  losses   increase  up  to 
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large  enough  values  to  offset  any  gains  realized  by  further 
increases  in  thruster  length.  These  viscous  loss  effects 
were  introduced  to  the  model  through  the  use  of  emperical 
data  used  in  the  formulation  of  the  expression  for  thruster 
voltage.  The  dominance  of  the  blowing  component  of  thrust 
can  also  be  observed  by  comparing  the  curves  for  thrusters 
10/5/1/. 006  and  10/2 . 5/ .  5/ . 0015  in  Figure  5.  These 
thrusters  are  the  same  length  but  the  10/2 . 5/ . 5/ . 0015 
thruster  one  half  the  radius  of  the  10/5/1/. 006  thruster. 
This  results,  again,  in  an  increased  average  current  density 
and  the  corresponding  decrease  in  efficiency  for  the  smaller 
diameter  thruster  for  constant  specific  impulse. 

Comparing  the  full  size  thrusters  to  their  scaled 
counterparts  discloses  the  scaled  thrusters  appear  more 
efficient  at  a  given  specific  impulse.  This  would  seem  to 
contradict  the  explanation  of  current  density  losses  in  the 
previous  paragraphs.  However,  examination  of  Figures  6  and 
7  show  that  the  efficiencies  of  subscale  thrusters  are 
quadratic  in  nature  as  opposed  to  more  linear  in  the  full 
size  thrusters  when  plotted  against  power.  A  continuation 
of  the  plots  for  the  scaled  thrusters  to  higher  powers 
results  in  crossing  of  the  full  size  and  subscale  thruster 
plots  making  subscale  thrusters  less  efficient  than  the  full 
sized  thrusters.  This  behavior  is  indicative  of  the 
presence  of  a  higher  order  loss  mechanism  beginning  to 
dominate  thruster  operation  at  high  specific  power  levels. 
This  higher  order  term  is  most  likely  a  reflection  of  frozen 
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flow  losses  being  introduced  with  the  empirical  elements 
used  in  the  formulation  of  the  computer  model.  Frozen  flow 
losses  are  the  result  of  energy  stored  in  the  excited  states 
of  the  propellant  as  it  is  exhausted.  As  the  propellant  is 
heated  and  ionized,  energy  is  stored  it  the  excited  states. 
The  propellant  is  exhausted  while  in  an  excited  state 
resulting  in  an  energy  loss.  The  density  of  MPD  thruster 
exhausts  are  too  low  to  have  collision  rates  high  enough  to 
de-energize  these  excited  states  and  recoupe  this  energy  in 
nozzles  of  any  practical  length.  Frozen  flow  losses  are 
related  to  the  energy  density  of  the  propellant.  Smaller 
thrusters  operating  at  approximately  the  same  power  have 
greater  amounts  of  energy  per  unit  propellant  mass  and  unit 
time  because  they  have  less  mass  in  which  to  deposit  the 
energy.  This  effect  is  present  throughout  the  range  of 
operating  conditions  but  is  more  pronounced  at  high  power 
due  to  the  increased  occurrence  of  secondary  ionizations. 

Figures  6  and  7  are  plots  of  efficiency  versus  input 
power  for  the  four  thrusters  examined.  These  graphs  show 
that  scaled  thrusters  can  be  expected  to  operate  at  lower 
power  levels  than  their  full  sized  counterparts  along  lines 
of  constant  efficiency.  This  phenomena,  as  the  above 
explanation  presents,  is  restricted  to  the  lower  power 
regimes  of  operation.  The  graphs  show  that  overall  power 
requirements  can  be  reduced  by  using  subscale  thrusters. 
What  is  also  shown  is  that  for  scaled  thrusters  operating 
along  lines   of  constant   efficiency,  the   smaller  thrusters 
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have  much  higher  specific  impulses.  This  is  not  necessarily 
an  additional  benefit.  Operationally,  a  propulsion  system 
has  a  minimum  total  mass  for  a  given  total  impulse 
requirement  determined  in  part  by  the  specific  impulse  [1]. 
Operating  at  specific  impulses  above  the  optimum  point  can 
result  in  an  increase  in  propulsion  system  mass.  This  is 
the  consequence  of  power  supply  mass  exceeding  the  optimum 
value.  The  graphs  show  a  one  half  reduction  in  linear  scale 
results  in  a  30%  to  40%  reduction  in  power  requirements  for 
thruster  operation  in  the  10%  to  20%  efficiency  regime. 
This  is  accompanied  by  a  40%  to  60%  increase  in  specific 
impulse  over  the  same  range  of  efficiencies. 
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IV.  CONCLUSIONS  AND  RECOMMENDATIONS  FOR  FUTURE  RESEARCH 

The  data  generated  by  the  computer  model  indicates  scaled 
thrusters  can  operate  at  comparable  efficiencies  to  full 
size  thrusters  but  at  reduced  power  levels.  The  reduction 
in  power  is  on  the  order  of  30%  to  40%  for  linear 
geometrical  scale  factors  of  one  half.  While  this  reduces 
the  power  requirements  for  an  operational  system  thereby 
making  it  more  plausible,  this  reduction  in  power  is 
accompanied  by  an  even  greater  increase  in  specific  impulse. 
The  increase  in  specific  impulse  is  on  the  order  of  40%  to 
60%  putting  the  subscale  MPD  thruster  into  the  2500  sec 
regime.  Studies  [8]  indicate  orbit  transfer  missions  are 
optimized  with  specific  impulses  in  the  1600  second  to  2000 
second  regime.  Scaling  of  MPD  thrusters  may  result  in 
operation  at  specific  impulses  higher  than  optimum,  thereby 
reducing  the  amount  of  payload  such  a  propulsion  system 
could  deliver.  This,  coupled  with  the  fact  that  this 
phenomena  only  occurs  in  the  undesirably  low  efficiency 
range  of  10%  to  20%,  indicates  scaling  thrusters  to  reduce 
power  requirements  may  not  be  a  good  choice  when  developing 
an  operational  system. 

Much  work  remains  to  be  completed  to  verify  this  model. 
This  will  require  extensive  testing  of  a  wide  variety  of 
different  MPD   thrusters  and   their  scaled  counterparts.   An 
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attempt  to  do  this  is  now  underway  in  the  Air  Force 
Astronautics  Laboratory's  Electric  Propulsion  Laboratory. 
However,  approximately  two  years  of  work  remain  before 
sufficient  empirical  data  is  generated  to  verify  and  further 
refine  and  expand  the  model.  Once  that  is  done,  the  model 
will  be  a  valuable  tool  for  the  prediction  of  thruster 
performance.  Appendix  A  is  a  data  reduction  program  written 
in  Fortan  and  intended  for  use  in  this  Laboratory  once  the 
test  apparatus  is  functional.  An  example  of  its  output 
follows  the  program. 
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n.±   ±  j-ili±jj.j\      .n 


DATA   REDUCTION   PROGRAM 


PROGRAM   MPD 

CHARACTER    HH*9,II*4 

REAL    V(40) ,U(40) ,W(40) ,VA(20) ,UA(20) ,WA(20) ,SI(40) ,SIA 
1(20),F(40),FA(20),E(40),EA(20) 

WRITE(*,106) 
106    FORMAT(lx, 'THIS    PROGRAM    ACCEPTS    UP    TO    33    DATA    POINTS') 
984    CONTINUE 

WRITE(*,101) 

101  F0RMAT(1X, 'ENTER    THE    CATHODE    LENGTH,    ANODE    LENGTH,     ' 
1,'(IN    CM),    AND    SCALE    FACTOR    IN    F4 . 1 ' ) 

READ(*,102)    XX.YY.ZZ 

102  F0RMAT(3F4.1) 
WRITE(*,001) 

001  F0RMAT(1X, 'ENTER    THE    DATE    AND    TIME    ACCORDING    TO    ' 
1,'AIR    FORCE    CONVENTION') 

READ(*,002)    HH,II 

002  F0RMAT(A9,A4) 
WRITE(*,981) 

981    F0RMAT(1X, 'ENTER    THE    MASS    FLOW    RATE    IN    KG/S    IN    F6.4') 

READ(*,086)    DM     - 
086    F0RMAT(F6.4) 

DMM=DM*1000. 

WRITE(*,199) 

199  F0RMAT(1X, 'HIT  CONTROL-P  THEN  RETURN  TO  PRINT  INTRO.' 
1, 'AFTER  PRINTING,  REPEAT  THE  COMMANDS  AND  REWIND  THE' 
2, '    PAPER    TO    THE   TOP. ' ) 

READ(*,198)    KP 
198    FORMAT(Al) 
WRITE(*,200) 

200  F0RMAT(9X, ' 
WRITE(*,201) 

201  F0RMAT(9X, ' 
WRITE(*,202) 

202  F0RMAT(9X, ' 
WRITE(*,203) 


203  F0RMAT(9X, ' | 
WRITE(*,204) 

204  F0RMAT(9X, ' | 
WRITE(*,003)    HH.II 

003    F0RMAT(9X, ' | 
WRITE (*, 205) 

205  F0RMAT(9X, ' | 
WRITE(*,206) 

206  F0RMAT(9X, ' | 
WRITE(*,207) 

207  F0RMAT(9X, ' | 
WRITE(*,208) 

208  F0RMAT(9X, ' | 
WRITE(*,209) 


THIS      DATA      WAS    OBTAINED    FROM 

THE        VARIABLE         GEOMETRY      MPD 

THRUSTER         AT      THE      AIR      FORCE 

ROCKET      PROPULSION      LABORATORY 

I 
ON       '  ,A9,  '       AT       '  ,A4,  '    HOURS. 

VOLTAGE    (V)         AND      CURRENT    (J) 

MEASUREMENTS    WERE    READ    FROM    AN 

OSCILLOSCOPE      ATTACHED      TO    THE 

THRUSTER    ELECTRODES.    THE    TOTAL 
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209  F0RMAT(9X, ' | 
WRITE(*,210) 

210  F0RMAT(9X, ' | 
WRITE(*,211) 

211  F0RMAT(9X, ' | 
WRITE(*,212) 

212  F0RMAT(9X, ' | 
WRITE(*,213) 

213  F0RMAT(9X, ' | 
WRITE(*,214) 

214  F0RMAT(9X, ' | 
WRITE(*,215) 

215  F0RMAT(9X, ' | 
WRITE(*,216) 

216  F0RMAT(9X, ' | 
WRITE(*,217) 

217  F0RMAT(9X,  '  | 
WRITE(*,218) 

218  F0RMAT(9X, ' | 
WRITE(*,219) 

219  F0RMAT(9X, ' | 
WRITE(*,220) 

220  F0RMAT(9X, ' | 
WRITE(*,221) 

221  F0RMAT(9X,  '  | 
WRITE(*,222) 

222  F0RMAT(9X,  '  | 
WRITE(*,223) 

223  F0RMAT(9X, '  | 
WRITE(*,224) 

224  F0RMAT(9X,  '  | 
WRITE(*,225) 

225  F0RMAT(9X,  '  | 
WRITE(*,226) 

226  F0RMAT(9X,  '  | 
WRITE(*,227) 

227  F0RMAT(9X,  '  | 
WRITE(*,228) 

228  F0RMAT(9X,  '  | 
WRITE(*,229) 

229  F0RMAT(9X,  '  | 
WRITE(*,230) 

230  F0RMAT(9X, ' | 
WRITE(*,231) 

231  F0RMAT(9X,  '  | 
WRITE(*,400) 

400  F0RMAT(9X,  '  | 
WRITE(*,234) 

234  F0RMAT(9X,  »  | 
WRITE(*,236) 


IMPULSE  (TI)   WAS  TAKEN  FROM  A 

TIME  INTEGRATED  ACCELEROMETER . 

TOTAL  GAS  PULSE  IMPULSE  (TIG) 
WAS  OBTAINED  FROM  FIRING  THE 
THRUSTER  WITH  THE  ELECTRODES 
SHORTED.  THE  DURATION  OF  THE 
THRUST  PULSE  (DT)  CAME  FROM 
MEASURING  THE  TIME  BETWEEN  THE 
HALF  CURRENT  POINTS  ON  A  PLOT 
OF  CURRENT  VS  TIME.  THE  RISE 
AND  FALL  OF  THE  CURRENT  WAS 
FAST  ENOUGH  TO  JUSTIFY  THE 
SQUARE  PULSE  SHAPE  ASSUMPTION 
INHERENT  IN  THIS  PROCEDURE. 
THE  MASS  FLOW  RATE(DM)  WAS 
DETERMINED  BY  MEARURING  THE 
SHAPE  OF  THE  GAS  PULSE  WITH 
AN  ION  GAUGE  AND  COMBINING 
THIS  WITH  THE  MASS  LOST  FROM 
A  PLENUM  TO  THE  GAS  VALVE. 

CALCULATED  DATA  WAS  OBTAINED 
FROM  THE  MEASURED  DATA  USING 
THE  FOLLOWING  EQUATIONS: 

THRUST(T)  =  (TI-TIG)/DT 
SPECIFIC  IMPULSE  =  T/(9.8*DM) 
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238 


239 


240 


241 


242 


243 
836 


009 


010 


Oil 


012 


013 


014 


015 


950 


951 


016 


WHICH  FITS   A 
XX,YY,ZZ,DMM 
THEM.  THRUSTER 


IK 


236  F0RMAT(9X, ' | 
WRITE(*,238) 
F0RMAT(9X, ' | 
WRITE(*,239) 
F0RMAT(9X, f | 
WRITE(*,240) 
F0RMAT(9X, ' | 
WRITE(*,241) 
F0RMAT(9X, ' | 
WRITE(*,242) 
FORMAT( '+ 

READ(*,243) 

FORMAT(Al) 

CONTINUE 

1  =  1 

L  =  0 

CONTINUE 

WRITE(*,010)  I 

F0RMAT(1X, 'ENTER 

READ(*,011)  V(I) 

F0RMAT(F5.1) 

WRITE(*,012)  I 

F0RMAT(1X, 'ENTER 

READ(*,013)  U(I) 

F0RMAT(F6.0) 

WRITE(*,014)  I 

F0RMAT(1X, 'ENTER 

READ(*,015)  W(I) 

F0RMAT(F5.3) 

WRITE(*,950) 

F0RMAT(1X, 'TO  CORRECT 
1 , 'ENTER  A  1  AND  ENTER 
2, 'ENTER  A  0.  ') 

READ(*,951)  JM 

FORMAT(Il) 

IF(JM.EQ.l)  GO  TO  009 

WRITE(*,016) 

F0RMAT(1X, 'IF  THIS  IS 
1 , 'GROUPING  OF  DATA 


EFFICIENCY  =  T**2/( 2*V*J*DM) 
AVERAGES   OF  REPEATED  FIRINGS 


ARE  SENT  TO  A  PLOTTING  PROGRAM 


SPLINED  CURVE  TO 


#'  ,3F4.1,F2.0, 


VOLTAGE  DATA  POINT  NUMBER  ',12) 


CURRENT  DATA  POINT  NUMBER  ',12) 


TOTAL  IMPULSE  DATA  POINT  NUMBER 


12) 


ANY 
THE 


OF  THE  LAST 
NEW  NUMBERS 


THREE  VALUES, 
OTHERWISE  ' 


THE 
POINTS 


LAST  OR  ONLY  ELEMENT  OF  A 
TO  BE  AVERAGED  TO  MAKE  ' 


2, 'A  SINGLE  POINT  FOR  PLOTTING 

3, 

4,'  ELEMENT 

READ(*,017) 
017  FORMAT(Il) 

IF(NN.EQ.O) 

L=L+1 

VA(L)=0.0 

UA(L)=0.0 

WA(L)=0.0 


'  ELEMENTS  IN  THE  GROUP.  IF  IT 
IN  A  GROUP  ENTER  A  0 
NN 

GO  TO  20 


ENTER  THE  NUMBER  OF' 
IS  NOT  THE  LAST' 


3^ 


M=I+1-NN 

DO  018  J=M,I 

VA(L)=VA(L)+V(J) 

UA(L)=UA(L)+U(J) 

WA(L)=WA(L)+W(J) 
018  CONTINUE 

VA(L)=VA(L)/NN 

UA(L)=UA(L)/NN 

WA(L)=WA(L)/NN 

WRITE(*,600) 
600  F0RMAT(1X, 'THE  AVERAGES  FOR  THIS  GROUP  ARE*) 

WRITE(*,450)  VA(L),UA(L),WA(L) 
450  FORMAT ( IX, Fl 0.5, 'VOLTS' ,F10.2, 'AMPS' ,F10.5, 'NEWTON-S' 

1, 'ECONDS') 
020  CONTINUE 

WRITE(*,30) 
030  F0RMAT(1X, 'ENTER  A  1  IF  YOU  HAVE  MORE  DATA  POINTS.' 
1, 'ENTER  A  ZERO  IF  YOU  ARE  DONE.') 

READ(*,040)  K 
040  FORMAT(Il) 

IF(K.EQ.O)  GO  TO  045 

1  =  1+1 

GO  TO  009 

045  CONTINUE 
MM=I 
WRITE(*,046)  MM,L 

046  F0RMAT(1X, 'YOU  HAVE  ENTERED  ',12,'  DATA  POINTS' 
1,'  CONSISTING  OF  ',12,'  GROUPS.') 

831  CONTINUE 

WRITE(*,050) 
050  F0RMAT(1X, 'ENTER  THE  TOTAL  IMPULSE  FOR  THE  GAS  PULSE') 

READ(*,060)  TIG 
060  F0RMAT(F5.4) 
WRITE(*,070) 
070  FORMAT( IX, 'ENTER  THE  TIME  BETWEEN  THE  ONE-HALF  * 
1, 'CURRENT  POINTS  FOR   A  TYPICAL  FIRING   OF  ' 
2, 'THETHRUSTER') 
READ(*,080)  DT 
080  F0RMAT(F6.4) 
WRITE(*,085) 
085  F0RMAT(1X, 'ENTER  A  CONTROL-P,  THEN  ENTER  A  BLANK  TO' 
1,'  PRINT  THE  DATA.  AFTER  THE  DATA  IS  PRINTED,  DO  THE' 
2,'  SAME  THING  TO  TURN  THE  PRINTER  OFF.') 
READ(*,089)  NP 
089  FORMAT(Al) 

WRITE(*,103) 
103  F0RMAT(46X,' ' 

1,' ~) 

WRITE(*,104)  XX,YY,ZZ 
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104  F0RMAT(45X,  *  | CATH0DE= ' ,F4.1, '  CM  |AN0DE=! 
1,  F4.1,'  CM | BENCHMARK  SCALE= ' ,F4 . 1 , * | ' ) 

WRITE(*,105) 

105  F0RMAT(*+ 

1,'      

2,' ') 

WRITE(*,350) 

350  F0RMAT(46X,' 

1  t  ' '  ) 

WRITE(*,351) 

351  F0RMAT(45X, ' |  GAS 
1 , ' | MASS  FLOW  RATE 

WRITE(*,352) 

352  FORMAT(  '  + 

1,*      

2,' 


IMPULSE 

D 


THRUST  PULSE  DURATION' 


35 


') 
WRITE(*,353)  TIG,DT,DM 


3  FORMAT (45X, '  I  ' , IX  ,  F5 . 4 , ' (N-S )  '  , 2X  ,  ' 
1, *(S) ' ,6X,  '  |  ' ,1X,F6.4, '(KG 


/S)\2X,» 


1 ,6X,F6.4, 
') 


WRITE(*,354) 
354  FORMAT( '+ 

1,'      

2/ 

WRITE(*,300) 
300  F0RMAT(46X, ' 
1/ " 


') 


.') 


WRITE(*,087) 
087  F0RMAT(45X  ' | VOLTAGE | AMPERES | TOTAL  IMP. (N-S) 
1,  '  ISP  |EFFICIENCY(%) |  ' ) 

WRITE(*,500) 
500  FORMAT( '+ 

1,'     _ 
2,' 

DO  090  1=1, MM 


.') 


F(I)=(W(I)-TIG)/DT 

E(I)=100.0*F(I)**2.0/(V(I)*U(I)*DM*2.0) 

SI(I)=F(I)/(DM*9.8) 

WRITE (*, 088)  V(I),U(I),W(I),SI(I),E(I) 


'  ,F5.1,1X, '  I  •  ,1X,F6.0, ' I  ' ,4X,F5.3,6X, 
5X,F4.1,4X, ' | ') 


088  F0RMAT(45X, ' 

l,'r,F5.0,' 
090  CONTINUE 

WRITE(*,301) 
301  FORMAT(,+ 

1,'     

2,' *) 

READ(*,093)  NA 
093  FORMAT(Al) 

DO  100  IN=1,L 

FA(IN)=(WA(IN)-TIG)/DT 

EA(IN)=100.0*FA(IN)**2.0/(VA(IN)*UA(IN)*DM*2.0) 


35 


100 


120 


130 


SIA(IN)=FA(IN)/(DM*9.8) 

CONTINUE 

WRITE(*,120) 

F0RMAT(1X, 'THE  AVERAGE  GROUP  VALUES  ARE.') 

DO  150  1=1, L 

WRITE(*,130)  VA(I),UA(I),WA(I),FA(I),SIA(I),EA(I) 


,1X,F6.0, * 
,4X,F4.1) 


,5X,F5.3,8X, 


1  IF 


YOU  NEED  TO  ENTER  A  NEW') 

') 
A  3  TO  GO  ON') 


831 
836 


F0RMAT(1X,F5.1,1X, ' 
1,F5.1,2X,' | \F5.0, * 
150  CONTINUE 

READ(*,700)  NB 
700  FORMAT(Al) 

DO  979  111=1, L 

UA(III)=UA(III)/1000. 

979  CONTINUE 
WRITE(*,832) 

832  F0RMAT(1X, 'ENTER  A 
WRITE(*,833) 

833  F0RMAT(1X, 'TIME  OR  MASS  FLOW  RATE 
WRITE(*,834) 

834  F0RMAT(1X, 'ENTER  A  2  TO  START  OVER 
READ(*,835)  NNN 

835  FORMAT(Il) 
IF(NNN.EQ.l)  GO  TO 
IF(NNN.EQ.2)  GO  TO 

980  CONTINUE 
WRITE(*,841) 

841  F0RMAT(1X, 'ENTER  990099 
WRITE(*,842) 

842  F0RMAT(1X, 'ENTER  965030  FOR 
WRITE(*,985) 

985  F0RMAT(1X, 'ENTER  960030  FOR  A  C0M1  PLOTTER') 
READ(*,843)  10, MO 

843  F0RMAT(I4,I2) 
IF(IO.EQ.9900)  10=10/100 
WRITE(*,972) 

972  F0RMAT(1X, 'ENTER  THE  LINE  SEGMENT  LENGTH') 
READ(*,973)  Z 

973  F0RMAT(F6.4) 
VA(L+1)=0.0 
VA(L+2)=60.0 
UA(L+1)=0.0 
UA(L+2)=6.0 
EA(L+1)=0.0 
EA(L+2)=10.0 
SIA(L+1)=0.0 
SIA(L+2)=1000.0 
CALL  PLOTS(0,IO,MO) 

CALL  WINDOW(0. ,0. ,10.7,7.5) 
CALL  PL0T(2.6,1.3,-3) 


2X 


FOR  CRT  DISPLAY') 

A  COM2  PLOTTER') 
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CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 


SYMBOL (-.4, 4. 3, .16,CHAR(0) ,90.0,-1) 
SYMBOL( -1.15 ,4.2, .16, CHAR( 2) ,90.0,-1) 
S YMBOL( - 1.9, 5.1, .16, CHAR (1), 90. 0,-1) 

=  ',90.0,3) 


,-13,7.0,0.0,UA(L+l) ,U 


986 


SYMB0L(-1.85,4.5,.16, 

LINE(UA,VA,L, 1,-1,0) 

CURVE(UA,VA,L,Z) 

LINE(UA,EA,L,1,-1,2) 

CURVE(UA,EA,L,Z) 

LINE(UA,SIA,L, 1,-1,1) 

CURVE(UA,SIA,L,Z) 

STAXIS(.l, .16, .1, .05,-1) 

AXIS(0. ,0. , ' CURRENT (KA) ' 
lA(L+2)) 

CALL  AXIS(0. ,0.  , 'VOLTAGE(VOLTS)  =  ! , 17 , 6 .0,90.0 , VA( 
lL+l),VA(L+2)) 

CALL  AXIS(-.75,0. , ' EFFICIENCY(%)  =  ' , 16 , 6 .0 , 90 .0 ,EA 
l(L+l),EA(L+2)) 

CALL  AXIS(-1. 5, 0. , 'SPECIFIC  IMPULSE ' , 16 , 6 . 0 , 90 . 0 , SIA( 
lL+l),SIA(L+2)) 
SIA(L+2)=500. 
SI(MM+1)=0.0 
SI(MM+2)=500. 
E(MM+1)=0.0 
E(MM+2)=10.0 
CALL  PLOT(0. 0,0. 0,-999) 
CALL  WINDOW(0. ,0. ,10.7,7.5) 
CALL  PL0T(1.3,1.4,-3) 

CALL  SYMBOL (-.4, 4. 3, . 16,CHAR(2) ,90.0,-1) 
CALL  LINE(SI,E,MM,1,-1,2) 
CALL  CURVE(SIA,EA,L,Z) 
CALL  STAXIS(.l, .16, .1, .05,-1) 
CALL  AXIS(0. ,0. , 'SPECIFIC  IMPULSE' 
CALL  AXIS(0. ,0. , 'EFFICIENCY(%)  = 
CALL  PLOT(0. 0,0. 0,999) 
WRITE(*,986) 
F0RMAT(1X, 'ENTER  0  IF  DONE,  1  IF  NOT') 


,-16,8. ,0. ,0. ,500.) 
,16,6. ,90. ,0. ,10.) 


READ(*,987) 
987  FORMAT(Il) 

IF(LLL.EQ.O) 

GO  TO  984 
983  CONTINUE 

STOP 

END 


LLL 


GO  TO  983 
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APPENDIX  B 
THEORETICAL  MODEL  PROGRAM 


001 
010 

005 

006 

007 

008 

020 

030 

040 

050 
056 

060 

061 
055 
070 


PROGRAM  MPDI 
REAL  J 

F(X)=1.00E-7*(J**4.0/X**2.0-4.94E-12*J**6.0*(X**2.0-2 
1.0*XMIN*X+XMIN**2.0))**0.5 

FF(X)=1.00E-7*(J**4.0/X**2.0-4.94E-12*J**6.0*(XMAX**2. 
10-2.0*XMAX*X+X**2.0))**0.5 
CONTINUE 
WRITE(*,010) 
F0RMAT(1X, 'THIS  PROGRAM  CALCULATES  THE  SUM  OF  THE' 


BLOWING,  PUMPING,  AND  ELECTROTHERMAL 


THRUST 


2, 'COMPONENTS  FOR  AN  MPD  THRUSTER') 
WRITE(*,005) 

F0RMAT(1X, 'ENTER  ANODE  LENGTH  IN  F4.1  CM') 
READ(*,006)  Z 
F0RMAT(F5.2) 
WRITE(*,007) 

F0RMAT(1X, 'ENTER  THE  MASS  FLOW  RATE  IN  KG/S  IN  F7.5') 
READ(*,008)  W 
F0RMAT(F8.6) 
WRITE(*,020) 

F0RMAT(1X, 'ENTER  THE  CATHODE  RADIUS  IN  F4.1  CM') 
READ(*,030)  XMIN 
F0RMAT(F5.2) 
WRITE(*,040) 

F0RMAT(1X, 'ENTER  THE  ANODE  RADIUS  IN  F4.1  CM') 
READ(*,050)  XMAX 
FORMAT(F5.2) 
F0RMAT(F7.1) 
WRITE(*,060) 

F0RMAT(1X, 'ENTER  AN  EVEN  NUMBER  OF  INTERVALS  IN  1100') 
READ(*,070)  N 
N=10 

CONTINUE 
WRITE(*,055) 

F0RMAT(1X, 'ENTER  AMPERES  IN  F7.1') 
READ(*,056)  J 
FORMAT(IIOO) 
H=(XMAX-XMIN)/(2.0*N) 
A  =  0.0 
AA  =  0.0 
AAA=0.0 
X=XMIN 

XXX=(XMAX+XMIN)/2.0 
DO  100  1=1, N 
XX=X+H 
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XXXX=XXX+H 

A=A+H*(F(X)+F(XX))/2.0 
AA=AA+H*(FF(XXX)+FF(XXXX))/2.0 
X  =  XX 

XXX=XXXX 
100  CONTINUE 
AAA=A+AA 
WRITE(*,130)  AAA 

130  F0RMAT(1X, 'BLOWING  THRUST  =  ',F8.4) 
A=XMAX 

B=XMIN 

Y=(XMIN+XMAX)/2.0 

AAAA=4.44E-13*(7.0*Y**3.0/6.0-(B**3.0+A**3.0)/3.0-(A*B 
l**2.0+B*A**2.0)/4.0-(Y**2.0/2.0)*(A*LOG(Y/A)+B*LOG(Y/B 
2))+(A*B**2.0*L0G(B/A)+B*A**2.0*L0G(A/B))/2.0)*J**3.0/Z 

WRITE(*,131)  AAAA 

131  F0RMAT(1X, 'PUMPING  THRUST  =',F8.4) 
AAAAA=3.25E-7*J**1.5*W**0.3*3.14159*(XMAX**2.0-XMIN**2 

1.0) 
WRITE(*,132)  AAAAA 

132  FORMAT( IX, 'ELECTROTHERMAL  THRUST  =',F8.4) 
AX=AAA+AAAA+AAAAA 

WRITE(*,133)  AX 

133  F0RMAT(1X, 'TOTAL  THRUST  =',F8.4) 

A1  =  (55.4-26.31*(1-XMAX/5.1')-23.4*(1-Z/21.6)) 
Bl  =  (39.214-65.84*(l-XMAX/5:.  1 )-27 . 656*( l-Z/2 1 . 6) ) *lE-4 
Cl=(35.7-54.5*(l-XMAX/5.1)+3.88*(l-Z/21.6))*lE-8 
V=(A1-B1*J+C1*J**2)*(178.-47.3*1000.*W+4.3*(1000.*W) 
l**2)/49. 
WRITE(*,200)  V,J 

200  F0RMAT(1X, 'VOLTAGE  =', F9 . 4 , 2X ,' CURRENT  =',F7.1) 
YX=AX/(W*9.8) 

WRITE(*,134)  YX 

134  F0RMAT(1X, 'SPECIFIC  IMPULSE  =',F8.3) 
E=( . 5*W*(YX*9.8)**2.0)/(V*J) 
WRITE(*,201)  E 

201  F0RMAT(1X, 'EFFICIENCY  =',F8.4) 
P  =  V*J 

WRITE(*,222)  P 
222  F0RMAT(1X, 'POWER  =',F9.1) 

WRITE(*,140) 
140  F0RMAT(1X, 'ENTER  1  TO  START  OVER,  2  TO  CHANGE  CURREN' 
l'T  ONLY') 

READ(*,150)  M 

150  FORMAT(Il) 
IF(M.EQ.l)  GO  TO  001 
IF(M.EQ.2)  GO  TO  061 

151  CONTINUE 
STOP 
END 
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